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O 

*j ' Abstract. We review Fujiwara's scheme, a sixth order weak approximation scheme for the 

i numerical approximation of SDEs, and embed it into a general method to construct weak ap- 

proximation schemes of order 2m for m 6 N. Those schemes cannot be seen as cubature schemes, 
but rather as universal ways how to extrapolate from a lower order weak approximation scheme, 
namely the Ninomiya-Victoir scheme, for higher orders. 

m 

1. Introduction 

Qh The Ninomiya-Victoir scheme for the weak approximation of solutions of stochastic differential 

Ph ' equations can be described in the following framework: let (Q, P) be a probability space 

and let . . . , Bf }t£R + be a d-dimensional standard Brownian motion. Define B® := t and 
Bt := (B® , Bl , . . . , Bf ). We consider stochastic differential equations driven by the Brownian 
motion {B t }teR+ 



C3 



d 



(1) X(t,x) = x + V [ V l (X(s,x))odB\ 

>. 

where x is in R N , Vi E C^°(R N ;R N ) and o stands for Stratonovich integral. We associate for 
later use the following simple stochastic differential equations to equation fl} 

(2) X^(t,x)=x+ I Vi(X^(s,x))odB, 



Let {Pt}teR + and {P t }ten. + be the associated heat semigroups on C£°(R d ) such that Ptf(x) := 

E[f(X(t,x))] for t > 0, and P t W f{x) := E[f{X^(t,x))} for t > 0. Notice here that the equation 
associated to the index is a pure drift equation, the semigroup a transport semigroup. Denote 
; furthermore by 

H ' i d 

»=i 

• - fp (0) o...oP (<1) 
Q? ] :=(P$o...oP$ 

Q? ] : = l(Qf ] + Q? ] )- 

the generator of the diffusion process (Q]), two ordered products of (semi-)flows with generators 
Vb and V? and the average of the two ordered products Q^. Then we have the well-known short 
time asymptotics, formulated in the language of fc-norms (see Definition [4]) 

\P t g(x) - Q? ] g(x)\ < Ct 3 \\g\\ 6e{d+1) , 

as t — » 0, leading - by iteration - to the Ninomiya-Victoir scheme. Indeed, when we define n-fold 
iteration of the operator Q ^ 

Q^„ = Q| ] o...oQg], 
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we obtain a scheme of weak approximation order r = 2, i.e., 

C 

\Prg(x) - Q [ T ] tn g(x)\ < ^\\g\\ 6 e(d+i)- 

Let us define formally weak approximations of Pt for a some fixed, finite T S R+ of weak 
approximation order r. 

Definition 1 (scheme of weak approximation order r). A family of linear operators {QT,n}ngN 
on C£°(R d ), continuous with respect to the supremum norm topology, is called a scheme of weak 
approximation order r if there exists C > and some number k > such that 

(3) \P T f{x)-Q T ,J{x)\<^\\f\\ k 

for all x € R N and for all / € C b °°(R d ). 

Notice that the operator Qt,u is only supposed to be linear and continuous with respect to the 
supremum norm topology on the set of C£°-function, but not necessarily of sub-Markovian type. 
This means in particular that classical (Romberg-)extrapolations belong to this class. 

In [5] T. Fujiwara constructs a sixth order scheme for smooth functions Cf°(R ) which consists 
of a linear combination of the previously described Ninomiya-Victoir scheme. Through the linear 
combination T. Fujiwara can "extrapolate" the weak approximation order to r — 6. In this paper, 
we define generalized Fujiwara schemes of order r = 2m including the scheme in [5] by refining 
Fujiwara's technique to prove the convergence order and construct versions of weak approximation 
order r = 2m for m G N. We finally obtain the following Theorem |4j whose proof can be found 
in Section^ notations can be found in the subsequent sections: 

Let {q n }n€ist be a generalized Fujiwara scheme of order 2m, then 

TCI 

Qt,u :=^/e,(Q| ] ) n 

i=l 

for n > is a scheme of weak approximation of order 2m, where a choice of k is given by 

m 

k = 2(2m + l)(d + l)^0i, 

i=i 

that means 

C 

\P T g{x) - Q T , n g{x)\ < -2^\\g\\k 

for test functions g € C%°(R N ). 

The remainder of the article is organized as follows: in Section 2 we introduce all algebraic 
prerequisities, in Section 3 we show the main algebraic result of this article, which is then applied 
in Section 4 to prove the existence of generalized Fujiwara schemes. In Section 5 we provide an 
implementation result, where the results can be compared to [10] . The appendix is devoted to an 
original proof of Fujiwara's basic algebraic result. 

2. Algebraic prerequisites and their relation to weak approximation 

Let A be a set whose elements are ao, • • ■ , as.- We call A an alphabet and do, . . . , letters. A 
word in alphabet A is a finite sequence of letters. Let 1 be a empty word and A* a set of words 
including 1. If we impose a total ordering on A, then A* together with word concatenation and 
lexicographic ordering becomes an ordered unital semigroup. Let R(A) be a set of noncommutative 
polynomials on A* over R i.e. a set of R-linear combinations of elements of A* and let R((A)) 
be a set of noncommutative series of elements of A* with coefficients in R, i.e. a set of functions 
/: A* — * R with well ordered support. Using componentwise addition and multiplication, which is 
induced by word concatenation, makes R((A*)) a R-algebra (see [4] for more details). The degree 
of a monomial is a number of letters contained in the monomial and the degree of a noncommutative 
polynomial and a noncommutative series are the maximum degree of monomials contained in them. 
Let H(A) m and R(A)< m be the set of homogeneous polynomials of the degree m and the set of 
polynomials of the degree less or equal to m respectively. Define H({A)) m and R((A))< m in the 
same manner. Since every u e R((^4)) has a well ordered support, we can define R((A)) >m ={a£ 
R((A))| deg(inf(supp(w))) > m} and R((A))> m = {u e R((A))\ deg(inf (supp(u))) > m} and it 
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is easy to see that R((A)) >m and R((A))> m are double sided ideals in algebra R((A)). Let j m 
and j< m be the natural surjective maps from R((A)) onto "R((A)) m and R((A))< m respectively. 

Since every subset of A* has a least element regarding lexicographical ordering, we have 
R((j4)) = H A . The set A* is countable, therefore taking metric topology in R makes H A 
with induced product topology into a Polish space. Hence, we can consider its Borel cr-algebra 
B(R((A))), R((A))-valued random variables and expectations, and other notions as usual. 

For u E R((A)) we define the exponential map 

it" 

exp («) ~22^> 

n>0 

and for u € H((A)) with vanishing constant term, we define the logarithm, 



Iog(l + tt) :=]T 



n 

n>l 



It is easy to check that 

(4) log (exp (u)) = u, 

(5) exp (log («)) = u, 

on the respective domains. For 8 e N define 



P : = exp(^aO, 



i=0 

exp(ia )---exp(-a d ) 

exp(ia d ) ■ ••exp(-a ) 
g «: = i(?i + Jfl). 

Let us make the substitution, which is the heart of the transfer from algebra to numerical 
schemes, a = Vq, a\ = . . . , a<j = formally correct. Let .B be another alphabet includ- 

ing vq, v\, . . . , Vd and set £?*, R(_B), . . . , in the same manner. For all t G R + define an algebra 
homomorphism ty t ■ R((.A)) — > R((_B)) by setting 

(6) * t (oo):=t«b, 

(7) %( ai ) := i«f/2. 

for all i € {1, . . . , d}. 

Define next an algebra homomorphism R(-B) — > C^ > °(R Ar ;R Ar ) by setting 

(8) *(«() = V5. 

Let -D = {X)«,eB* a ™H Stoes* a ™ <&(«>) i s wen defined }. Clearly, R(-B) C I? and Z? is a R- 
subalgebra of R( (£?)). The homomorphism <f> can then be uniquely extended to an R-algebra 
homomorphism D -► C b 0O (R JV ; R"). 

The algebra of non-commutative words plays a major role in the analysis of weak approxi- 
mation schemes due to the following well-known asymptotic expansion theorem, which allows to 
approximate the truncated exponential series in A by other simpler expressions. 

Theorem 1. For all function f <E C£°(R N ), x € R w and neN, it holds that 

n ,j. 

(9) PJ(x) = £ ^7(z) + 0(t n+1 ) = <f>(Mj<nP))f(x) + 0(t n+1 ). 

k=0 

as t —>■ 0. 

Proof. See [B]. □ 
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Hence we can, e.g., express the generator A of the diffusion process (fT|) by 

$(*i(oo + ... + o d )) =A, 
in particular we obtain the following crucial asymptotic formulas, 

$*tU<n(exp(ai))) = Pt i] + 0(t n+1 ) 

as t — > and i = 0, . . . , d again due to Theorem [TJ 

To be more precise on the goal of our paper, Theorem [T] also means that if we approximated p 
by linear combinations of (q^) up to a certain degree 2m — 1 within the algebra R(A) such that 
the remainder term is of order 0(-5sr), then P t f(x) could be approximated by linear combinations 
of $(^> f ((q^) n ))f(x) in a weak sense of order 2m. 

Notice that the letters at correspond to squares of vector fields under $ o hence one has to 
work out the correspondence to exponentials of first order terms, too. The next lemma shows how 
to relate thoes linear semi-flows of PDEs pj 1 ' to non-linear flows of ODEs Fl^ (x) up to a certain 
degree m, namely by replacing the normal random variable Z by a random variable taking finitely 
many values and sharing moments up to order 2m. This finally means that we can approximate 
by convex combinations of exponentials of first degree terms, i.e. ao, . .. ,dd leading to weak 
approximation schemes. 

Lemma 1. For all i £ {1, . . . , d} we have that 

(10) £[exp =exp(i|-) 

holds true. This formula also holds true under the homomorphism $ o v^, i.e., 

E[f(Fl^(x))]=P^f(x) 

for test functions f and x £ R w . 

Proof. Proof by applying the Fourier transform of Brownian motion and classical subordination 
results. □ 



3. HOW TO APPROXIMATE p BY q? 

An alternative proof of this result can be found in the appendix: 
Lemma 2 ( [5] Lemma 2.1). We have 

. OO y 

(11) log?™ ^(-ly+^loggW). 

i=l 

Proposition 1 ([5] Proposition 2.2). There exists Ci £ R((A))>2i+i such that for all 9 £ N, 

OO 

holds. 

Corollary 1. Let q be a linear combination of for some 9 £ N. If there exists n £ N such 

that j< 2n -l(q) = 3<2n-l(p)> then 3<2n{<l) = J<2n(p)- 

Proof. For all £ N, j<2(q' 6 '') = j<2(p) holds. Hence, the case n = 1 is clear. Suppose n > 2 and 
j<2n-i(?) = j<2n-i(p)- Since g = X)j=i "i?^' 1 for some e N and sincc j<2(<? [el ) = j<2(p) for 
all £ N, it follows 53j=i a 3 ' = 1- According to Proposition [T] 

V /C ^ 

i=l J=l J 

for some Cj € R((A))>2i+i. Since j'<2n-i(?) = i<2n-i(p)j we nave 

fc 1 

= ° 
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for all i — 1 , . . . , n — 1 . Then 

oo k 
i—n j—1 3 

which proves the corollary. 
Set 



A := 



1 

1/91 



1 

Vol 



1/6 



2(m-l) 



Corollary 2. 



(12) 



3<2r, 



( 


/ 


"1" 




T 




\ 











~ q Wi] - p - 






A- 1 












V 


V 


_0_ 


J 




_ g [fln.]'_p_ 


/ 



holds. 

Corollary 3. for alHg{l,...,m- 1}, 



(13) 



/ 


"1" 






T 


' 1 " 

ef 


A- 1 








i 


V 





J 




_ m _ 



□ 



4. Generalized Fujiwara scheme and its property 
Definition 2 (Generalized Fujiwara scheme). A family of series, 

m 

4 = 1 

is called a generalized Fujiwara scheme of order 2m if 

/ = [/ft ' ' ' fe m ] T 




holds. 

A straightforward calculation involving induction gives the following connection concerning the 
powers of series in R((^4)). Notice that we split the product q n — p n into telescoping summands, 
where one, two up to m terms of the form q — p appear. 
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Proposition 2. For p,q G R((A)) and 2 < m < n, we have 



k=0 

m — 1 n—l hi — I k-s — 1&2— 1 

+ E( E E ■•■EE? fcl ^P)P (w, (?-?) 

1=2 k,=l~l fe,_ 1 =i-2 fc 2 = lfci=0 

x p(*l-*»-l-l)(g_ p )p»-fcl-l) + 
n—l fc m — 1 fes— 1^2— 1 

E E -E E^-^-rfx- 

fe m =?7l— 1 fc^ — i— m — 2 ^2 — 1 fci=0 

~km—km — l — 1 ^n — fc T7X — 1 



x jr 



In particular for m = 1, 



l (q-p)p n 
n-1 

q n -p n = J2<l k (<l-P)P n ~ k ~ 1 

k=0 



holds true. 

Lemma 3. For Zi, z 2 G R((A)), if j<i(zi) = and j< m (z 2 ) = 0, t/ien j< i+m+ i(ziz 2 ) = 0. 

Proof. By the assumption, monomials with the lowest degree contained in z\ and z 2 are of the 
degree I + 1 and m + 1. Then, monomial with the lowest degree contained in ziz 2 has the degree 
I + m + 2. Hence j<> +rn+ i(.2iZ2) = 0. □ 

Corollary 4. Forzi,z 2 ,z 3 € R-((^)), if 3<l{z\) = i<i{z2) and j< m {z 3 ) = 0, then j<i +m +i(ziz 3 ) = 

j<l+ m + l(z2Z 3 ). 

Corollary 5. For z G R((A)), if j<i{z) = 0, then j< m/+m _i(z m ) = 0. 
Theorem 2. If a series, 

m 

(15) qn:=Y,^(q [e ' ] ) n , 

i=i 

is a 2m-th order generalized Fujiwara scheme, then for all I G {2, . . . , m — 1} , 

m 

(16) j< 2m+H (E/»,(9 |e ' L p) , ) = 

i=l 

holds true. 

Proof. Fix I G {2, . . . , m — 1}. By Proposition [U 

oo 

Czi ' ' ' 



d7) Efcfe 1 * 1 -*)' = £/* E 

z=l z=l ii,...,ij=l "z 



holds. It is easy to see that 
Hence, we have 



Cjj ■■■C il G R((A))>2( 41 + . ...,,)+;. 

(00 
\ ^ c ii c ii 
2^ fl 2(i 1 + ...+i,) 
H,...,i ! =l t'z 



=J<2m+/-l 



-J<2m+i-l 



E 



zi,...,ii>l 
\z'iH Ki<m— 1 



02(ii+.-i,) 



/ 



7/4 — ± 

E V c zi " " " 
Q2k 



k=l Zi,...,Z;>l 
z'l + -..+i;=fc 
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Thus, we have 



3<2m + i-i(^feM ei] -P) 1 ) 



1=1 

m— 1 



=J<2m+l- 



^E E 12fa-^k c ^--- c ^ 



k—l i\,...,ii>l i—1 

ilH Hi=fc 



m— 1 



k—l ii,...,i/>l 
iiH Hz=fc 



Or' 1 ' 



1 



=0 



by Corollary [3l 



□ 



Definition 3. A generalized power series a £ H((A}) is an element of O(s) if for every n, N £ N, 
n < N there exists a uniform bound for all the coefficients of the terms of -a with degree k, which 
satisfies rt<A;<iVass^0. 

Theorem 3. //{(j n }n€N is an m-th order generalized Fujiwara scheme, then 

m n— 1 fe^— 1 ^3 — 1/22—1 

*X/nqn-*lP = 12 12 12 '•• £ X] *l/n°i,(*i,....fci)i 
J=l fc(=i-l fe i _ 1 =i-2 fc 2 =lfci=0 

w/iere j<2m+i-i(a;,(fc i ,...,fc 1 )) = /or I = 1, . . . , to, and ^i/ n q n ~ £ O(^asr). 

Proof. Let m > 2. The case m = 1 is trivial. Let {g n := ^2=1 fo i (q^ e ^) n }ne'N be an m-th order 
generalized Fujiwara scheme. Note that ^>xP = {^\/ n p) n ■ Then by Proposition [2l we have, 

m 7i— X 

= 12 fa 12^^p) k ^l/nq [e ' ] - ^l/nP)^l/nP) n - k - 1 

i=i k=o 

m m — 1 n— 1 fc; — 1 £3— 1 &2 — 1 

+ E^E E E • • • E E (*l/nP) kl (*l/n<i m Vl/nP^ynP)^- 1 X • • • 
i=l i=2 fcj=i— 1 fe ( _i=i— 2 fc 2 = lfc 1= 

x (* 1/n9 ^] - *i /n p)(*i /n p) fc, - fe! - 1 - 1 (*i/„g M - * 1/nP )(* 1/ „p)"- fc '- 1 



n-l fcm— 1 



fea-l fc2-l 



+12 fa E E -EEfW^-VKV) 



^2 — fci — 1 



i—1 k m —rn — l k rn — \=m—2 k%=\ k\ — 

'x/nQ. 1 '" ~ 'i' i , 



X (* 1/n q m - *l/„p)(^l/nP) fcm_fcm - 1 ^l/n^* 1 - *X/nP)(*l/nPY 



Set 



ai, (fcl) = E^ fcl (« Ie<] -P>' 



n—ki — 1 



For / e {2, . . . , to - 1} set 



ai,(k l ,...,k 1 )=12faP kl (i m -p)P' 



k 2 -k 1 -X ( W i ] 



P)---P 



fcl— fcl-i — 1/ 



9 4 ]_ p ) p n-fc.-l 



and for / = to define 



am,(fc ra ,...,*0 = £/*(? [9il )* 1 (9 [9,] -p)** 9 "* 1 " 1 ^ 1 -P) ■ ■ -P 

i-X 

In particular the summand dxrki) can be written as 

m 

(18) a Ukl ) =P kl 12fa(<l le ' ] -p)P 



i=X 



km — km-X — Xl 



P)P' 



i-k m -X 



n—ki—1 
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Let Cj S R((A))>2j+i, i € N be as in Proposition [T] By Theorem 

m 

i<^(X>«(9 lfl ' ] -p)) = 0, 
»=i 

thus J<2m( a i,(fei)) = 0. Also it holds that 

rn 



By Corollary H for all 9 € N, j< 2 (<7 [01 - p) = holds. Thus, by Corollary EJ Corollary 0] and 
Proposition [1] 

(19) j< 3m -l((q m ) k H<i m - p^-^-Hl™ -P) ■ ■ V 
holds, hence j<?, m -i{a m ,(k m ,...M)) = °- Moreover, 



p)p»-*»-i) = 



Wfa 1 * 1 )* 1 ^ 9 * 1 -p)^-* 1 - 1 ^ 8 *] -p)---p fc " 



-km-l-ll 



p)p«-^-i) = ( j3 ( C i)) m ^r 



Let pi, . . . ,p;+i G R((A)) with property Jo(ft) = 1 for * € {1, . . . , Z + 1}. By using similar 
arguments as in the proof of Theorem [2 we get 

m 

j<2m+I-l(£/fl«Pl(« [ ° <] -P^i^ -P)---Pl(i m -P)PI+1 

PlCiiP2 ■ ■ -piCitPl+lj 



2=1 



3<2m+l-l 



EE/ 5 

k—l ii,...,i|>l 
iiH Hi = fc 



"r' 



0. 



and 



j2m+l(j2MQ [ei] ~P)P2(q m -P)---Pl(q [(> ' ] ~ p)P(l+l)) 



= E f 1 

ii,...,ij>l 
iiH Mi=m 



for alU € {2, . . . ,m - 1}. We conclude, that j<2m+l-i(ai,(k,,...M)) = °- 
It remains to prove that ^>i/ n q n — *ip £ 0{ n im+i )■ 

First, let us observe auki,...,ki) f° r ^ G {1, . . . ,m — 1}. Choose Me NU {0}. As above we can 
write 



j2m+l+M(0'l,(ki,...,k 1 )) 



-J2m+1+M 



Y,^p kl ^ [e ' ] -p)p k2 ~ kl ~ 1 ^ ] -p) 



■ ■ -p 



ki-ki-i-ltq[6i] _ p\pn-h-l 



m+[M/2] 

Y, E f 1 

k—m ii,...,i/>l 

— +*!=*! 



0? fc 



1 



J2m+l+M\P C tl p ■■■p C,,p J 



The coefficient of the power p k of the term of degree I is of the form k l c, where c is the coefficient 
of p of the same degree. Hence, the coefficient of the term of degree 2m + I + M of 



is a finite sum, namely 



r,kl „fc2— fcl-l fcj— — 1 71 — fc; — 1 



£ fci 11 " (*a - kx - l)" 2 - 1 . . . (n - fc, - l)" ! + 1 - i 6 i 
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where 6j and the number of summands do not depend on n, and rij^ € N, Y^jtX n j,i — ^ m + 
M — 2J2k=i Let us denote b' i k — . . . ^] fbi. Thus, the coefhcicnt of a term of degree 
2m + 1 + M of 

n— 1 fe;— 1 fe 3 — lfc 2 — 1 

E E ' ' ' E E a i,(fci,-, fci) 
ki=l-l ki-!=l-2 fe 2 = lfc 1= 

has the following upper bound 



^ n— 1 fc; — 1 fe3 — 1 fc 2 — 1 

,2m+M+i E E -EE 



ki=l-l fci_i=i-2 fe 2 = lfei=0 
m+[M/2] 

E E E fc " M ( fc 2 - fci - l)™ 2 '* ■ • • (n - fe, - 

fc=m ii,...,ii>l i£-T 
iiH Hi=k 

^ n — 1 fe;— 1 &3 — lfc 2 — 1 

E E --EE 

fe;=;-l kl-!=l-2 fe 2 =lfei=0 

m+[M/2] 

E E E^" 1 " ( fc 2 - fei - ir<* ...{n-h- ir+^b'j 

k—m i i , > 1 z£7 

iiH Hi=k 

n-1 fei-1 fe 3 -l fe 2 -l m+[M/2] _^ 

^E E -EE E E E ^euJ ^'- 

k t =l-l ki-!=l-2 fc 2 = l fei=0 fc=m ii,...,i,>l ie/" 

iiH Mi=k 

Since the number of terms of 

n— 1 — 1 fc3 — 1/C2 — 1 

i2m+;+M(*l/n E E '"EE a i,(fc|,-,fci)) 
fc;=i-l fe ; _i=i-2 fe 2 = lfci=0 

is finite and their number does not depend on n, there exists a uniform bound for all of the 
coefficients, which proves our assertion for I < m. 
Let us observe the a m ,(k m ,---M)- ^ s aDOve 

j2m+l+M(a m ,(k m ,. ..,fei)) 
m 

=jW J+M (E^(9 [fl4l ) fcl (9 [fl4] -P)p fa - fcl - 1 («^ -p)-/™-^'- 1 ^^-^- 1 ) 
i=i 

m+[M/2] m 

= E E E^^+Mi^i)^^/-^- 1 .../-^-- 1 ^"-^- 1 ) 

fc=m ii,...,i m >l i=l 4 

Since is a convex combination of products of exp(^aj) all its coefficients are positive and the 
sum of all coefficients at the terms which are derived from a^ ■ ■ ■ a s d d by permutation of dj's is 
exactly the coefficient at the term a^ ■ ■ ■ a s d d of the exp(J^ i=0 a-i) in the (commutative) power series 
algebra, generated by oo, ■ • ■ , ad, with coefficients in R and the same goes for the coefficients of 
the power (q^) k and the coefficients of the commutative series (exp(^ i=0 a,j)) k . The rest of the 
argument goes as in the case of I < m, which gives us the upper bound of the coefficients at the 
terms of 

n— 1 k m — l &3 — 1 fc 2 — 1 

j2m+l+M 

(*Vn E E E E a m,(k m ,...M))- 

km—m—1 k m _i—m — 2 k 2 —l fei = 

□ 



Now we are able - by means of our homomorphisms W and $ to transfer the algebraic results 
into the realm of weak approximation schemes. 



in 
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Definition 4. Let ieN and let g G C b °°(R d ). Define 

IMU := sup || v l fflloo- 

Remark 1. The function || • ||fe is a norm on C^°(R d ). 

Definition 5. Let denote the space of bounded linear operators on (C£°(R d ), || • ||fc). We can 
regard Bk as a normed space with the operator norm. 

Proposition 3 ([6]). Fix a k £ N. The following assertions hold: 

(1) The family (Pt)t>o is a uniformly bounded subset of Bk- 

(2) Let A be the geneartor of the diffusion process |T]) and let N G N. Then, for g £ C£°(R d ) 
we have 

N f k If* 
(PtffX*) = £ + j- / (t - ^(P.^sX*)?.^- 

We are now able to formulate weak approximation schemes of order 2m by means of our 
algebraic preparations. Recall therefore the definitions 

qW . _ f P (0) o .„ oP M 

Qf ] : = i(0F + SF)- 
of the building blocks of Ninomiya-Victoir schemes. 

Theorem 4. Let {cfaJneN fre a generalized Fujiwara scheme of order 2m, then 



Qr,n ■— ^/^(Qt 1 )" 

i=l 

/or n > is a scheme of weak approximation of order 2m, where a choice of k is given by 

m 

k = 2(2m + 1)^6^, 

i=l 

that means 



C 

\P T g{x) - Q T , n g(x)\ < -j^HsflU 



for test functions g £ C%°(R N ) 



Proof. Due to asymptotic formulas 

^t{]<2m{p))=Pt+0{t 2m+1 ) 

and 

** t (j< 2m (exp(a l ))) - P t W + 0(t 2m+1 ) 
where the constants in the Landau symbol depend on the derivatives of order at most 2 (2m + 1). 
Therefore we can simply copy the proof of Theorem [3] by first replacing q n with Qt,u and p by 
Pt- In the appearing sums we have to use the previous asymptotic formulas, namely 

+ ^I.U<2m((q le] ) n ~p)) + 
+ Wz(i<2m(p))--P£, 

where the order behavior of the middle part has been shown in Theorem [3] and the order behavior 
of the other two summands follows from the previous asymptotic formulas. Apparently each term 
in Q^ s \ which is approximated due to the asymptotic formulas, increases the number of derivatives 
necessary to do the estimation by 2(2m + 1), which leads to the formula for k. □ 

Example 1 . The case m = 1 apparently corresponds to a version of the original Ninomiya-Victoir 
scheme. 
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Example 2. The case m — 2 corresponds to a scheme already presented in [5]. One can choose 
6% = 1 and 6*2 = 2 and = — ^ and /e 2 = f ■ 

Example 3. The case m = 3 corresponds to Fujiwara's originally presented scheme, which in our 
language reads like follows. Notice that we do not need the full strength of our previous proof, 
which is built on Theorem [2] 

For all mutually different numbers 61,62,83 £ N, we can construct 6-th order generalized 
Fujiwara scheme q with a form: 

q = feM 6l] T + /fafe 1 * 1 )" + fez(q [93] ) n - 
For the proof, which is presented for convenience here, we assume without loss of generality 
that 61 < 6 2 < 63. We have, 



fe 2 
fe 3 



By Corollary [2] we have 



{Bi-ffiWi-ffi) 

(ei-ffiWi-ffi). 



K4(^- P )=Ja(|(? N 

M<P* - p) = j<«(§(ff [ * ] 
63 

Then, by Corollary 21 we have, 



62 

1 



p)), 
P))- 

-P)% 
-Pfl 



Thus, 



j<7(J2M(i m ~p) 2 ) 



i=l 



= ( 



(0 2 2 - 6()(6i-e() {9i~9i){6. 



= 0. 



(oi-oi)W-oi) 



)3<7((q l6l] - pf) 



5. IMPLEMENTATION OF A TO-TH ORDER GENERALIZED FUJIWARA SCHEME 

A scheme for approximation of expectation of order six was first introduced by Fujiwara [5]. In 
previous sections we theoretically constructed schemes for approximation of expectation of order 
2m for arbitrary m € N. In this section we show how to construct a practical scheme with 
approximating flow of vector fileds Vt, which drive the SDE |T]), by some suitable integration 
schemes. The usual choice for the integration schemes are Runge-Kutta methods. In our concrete 
example from mathematical finance we will use a seventh-order nine-stage explicit Runge-Kutta 
method with a very good stability, given by M.Tanaka et al. (see [11], [12] and [E]). Higer 
order Runge-Kutta mehod often lose stability with respect to rounding error, truncated error and 
piling error. In addition, these effect decrease order of approximating error. Since in a concrete 
application of the algorithm, e.g. in mathematical finance, some of the ODEs can be very close to 
being stiff, the stability of the Runge-Kutta algorithm is of high importance. We show a relation 
between convergence order of weak approximation scheme and m-th order Runge-Kutta method. 
In addition we construct a concrete algorithm of a m-th order generalized Fujiwara scheme and 
analyze its computational cost and its approximating error. At the end we present a concrete 
numerical experiment. Tanaka's result is presented in the Appendix since we could not find any 
of his papers written in English. 

The results of this section can be compared to those from [TO] . 
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5.1. Runge-Kutta method. For V E C£°(R N : R N ), the map exp : C£°(R N ,R N ) x R + x R N 
represents the flow driven by the vector field V starting at xq, i.e. the solution of the ordinary 
differential equation: 

(20) |*m=v(* (t )), 

x(0) = Xq. 

Definition 6 (s stage explicit Runge-Kutta method of order m for autonomous systems). A s 
stage explicit Runge-Kutta method of order m for autonomous systems is determined by a lower 
triangular matrix A = [a^]^ - =1 and a row b = [bi ■ ■ ■ b s ] such that the following hold: 

• Let h 6 R, to E R and let t n = t n -\ + h for all n 6 N. Given the vector x n -i as an 
approximation to x(i„_i), where x satisfies the equation (|20p. the approximation x n to 
x(t n ) is computed by evaluating, for i = 1, 2, . . . , s, 

Fi = V(Xi), 

where X\ , X 2 , . . . , X s are given by 

Xi = x n -i + h aij-Fj 



and then evaluating 



3<i 



Un = Vn-i + hJ~]bjFj. 

i=i 

• The Taylor expansion of function of h around should coincide with the Taylor 

expansion of x(t n ) = x(t n -i + h) up to (including) the term at the power h m+ . 

Remark 2. Usually Runge-Kutta methods are studied for general non-autonomous systems. In 
these cases the method is uniquely identified by a triplet A, b and c, where A and b are as above 
and c = [ci . . . c s ] T is a suitable column vector. 

See Butcher [2] and [5] for more details about the theory of Runge-Kutta method. 
The next theorem shows that we need at least 12-th order Runge-Kutta method for 3-rd order 
generalized Fujiwara scheme. 

Theorem 5. For all f E C%°(R N R N ) ; t E R+ and x E R N , there exists C t > such that 

\f(exp(tV )) - f(R m (t,V )(x))\ < C t m+ \ 
\E[f(exp((VtZVi)) - f(R 2m (VtZ,V l )(x))}\ < Qt m+1 , 
where i E {1, . . . , d} and Z ~ Af(0, 1). 

Proof. The first inequality follows from the definition of m-th order Runge-Kutta method and 
Taylor's theorem. Set i E {1, ...,d}. By the definition of Runge-Kutta method and Taylor's 
theorem again, we have, 

/(exp ((VtZV t )) - f(R 2m (V~tZV l )(x)) 

y-rn+l/2 ^2m+l 
= : - y2m+l if \ , Q(t m + 1 \ 

2(m+l)! K JW + u\? )■ 

Note that for all k E N, E[Z 2k+1 ] = holds. Thus the conclusion is true. □ 

The next theorem shows that if we do not urge to have O(n) computational cost, 4th order 
Runge-Kutta method is enough for sixth order scheme. 

Theorem 6. For k, n E N, for all f E C£°(R N ), for all i E {1, . . . d}, and for all x E R N , there 
exists Ci > such that 

\E[f fexp (-^vX) - f (R m (^,vX (x))}\ < ( ' 



km+k+m/2+1 



holds where Z ~ A/"(0, 1) 
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Proof. 



E\f(exp (-^=Vi)(x)) - /(i? m (-A=l, Vl f (x))}\ 

= \E[(e X p-^V i r k f(x)-R m (^,V i ) nk f(x)}\ 
n K \Jn n K \Jn 

n fc -l 

1=0 v v v 

R m {-^,v i r h - l - x f{x)] 

. Ci u 

< 77 

- n 2(fe+l/2)(m/2+l) 

< 

— n (km+k+m/ 2+1) ' 



a 



5.2. Recipe for to— th order generalized Fujiwara scheme. In the following subsection we 
will provide the pseudocode for implementation of the m-th order generalized Fujiwara scheme 
with fixed coefficients 9i, $2, ■ ■ ■ , m . Let / = [/i • • • f m ] T be as in the section [3] and let the function 
solveDE(V,xo,t) return the solution of the ODE (|20p at time t with initial condition x(0) = xq. 



Algorithm 1: Fujiwara 
Data: function g, vector fields Vo, V±, . . . , Vd, time T, initial condition xq, number of 

partition points n, number of samples M 
Result: approximation E of the expectation Fi[f(Xx)], where X t is a process defined by the 

SDE Q 
Q <- S R lxm ; 

for o <— 1 to M do /* expectation (MonteCarlo or quasi Monte Carlo) */ 

I Q <— Q + samplePath{g, Vq, . . . , Vd, T, xq, n); 
end 

Q < jfQ] 

/* approx. for E{g(X{T, Xq))) is the linear combination fi* Qi */ 

/••• Qf- 

return E 
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Algorithm 2: samplePath 



Data: function g, vector fields Vq,Vi, . . . , Vd, time T, initial condition xq, number of 
partition points n 

Result: row vector Q = [Qs ■ ■ ■ qI"" 1 '] G R lxm calculated for a random simulated path 
()■ () •: Ti '": 

generate independent Bernoulli(l/2) random variables A = [Ai, . . . , A n ]; 
for k <— 1 to to do 

/* generate all standard normal random variables that are needed */ 

D <— T /(nflfc)[l ■ ■ ■ 1] € R lx (^"); 

N <- y/T/(n0 k )[Ni, Ne kn \, where N x , . . . , Ng k7l are i.i.d, with Ni ~ iV(0, 
/* first row serves for the component without Brownian motion */ 
~D] 



[N 

X <- x Q ; 

for j <— 1 to n do 



if A,- = 1 then 



/* consequtively solve the DDE's */ 
/* solve appropriate ODE */ 
for <— 1 to 9k do /* repetition because of finer dissection */ 

for i <— to d do /* solving DDE's */ 

[ X «- solveDE(Vi, X, Z l+1:ie _ 1)m+J ); 
end 
end 
else 

for 9 <- 1 to fc do 



/* repetition because of finer dissection */ 



for i -e 

I X 
end 

end 

end 

end 



to d do 

- solveDE(Vd-i,X, ^d- 4 +i,(e-i)*«+i); 



/* solving DDE's */ 



end 

return Q 



Remark 3. Usually in modern computers memory size is no longer an issue. From this perspective 
it seems sensible to generate all needed random variables in advance. Namely, the random variables 
for various O^s do not have to be independent, therefore we can reduce its number by reusing them, 
and there exist efficient algorithms which speed up the process of their generation if we do it in 
one batch instead of step by step as it is written in Algorithm [2] 

5.3. Computational cost. 

Theorem 7. Let d,n, M,m,T,9i, . . . ,9 m be as above, such that T/n is sufficiently small. Fur- 
thermore, assume that each step of the method solveDE needs a operations, i.e. additions, multi- 
plications and function evaluations, that B operations are needed to generate a (pseudo or quasi) 
Bernoulli random variable and that Z operations are needed to generate a standard d-dimensional 
normally distributed (pseudo or quasi) random variable. Then the computational cost of Algorithm 

□ is MUm + n((d + l)o + Z + l) £™ =1 9 k + nB + lj+ 2m. 

Proof. Let us denote the computational cost of the Algorithm [2] by C. A straightforward calcula- 
tion shows that the computational cost of the Algorithm [T] is equal to M(C + 1) + 2to. 

For fixed j 6 {l,...,n} in Algorithm [2] we have 9 k {d + l)a operations. Hence, for fixed 
k E {1, . . . , to} there are 5 + n9k((d + l)a + Z + 1) operations. It follows that C — 5m + n((d + 
l)a + Z + 1) £fcL x °k + nB. □ 

Remark 4. Rigorous use of Runge-Kutta algorithms for solving ODEs in the algorithm is only 
suitable for building a universal solver for SDEs of the type (QJ . In concrete practical applications 
it is to be expected that many of the ODEs of the type (|20|) have a nice enough explicit solution 
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The error of the algorithm consists of discretization part, i.e. the error due to numerical solution 
of ODEs and the error which comes from the scheme, and of the convergence error which comes 
from the Monte Carlo or quasi Monte Carlo simulation. 

Theorem 8. For n, M S N such that T/n is sufficiently small. The approximation error of 
Algorithm 1 is 0(l/n 2m ) + 0(1/ VM). 

Remark 5. One should take great care when choosing a suitable subdivision of the interval, since 
the coefficient of the discretisation error directly depends on function / and vector fields V%, thus, 
although bounded, the coefficient can get fairly large in some cases. Moreover, the convergence 
error of the Monte Carlo simulation is directly proportional to the sqare root of variance of 
f(X(T, x)). As in the case of discretisation error this should be taken into account, since, although 
constant, the variance can be large comparing to the size of error we would like to achieve. 

5.4. Numerical example. For our numerical example we have chosen the genearlized Fujiwara 
scheme of order 6 with 6\ = 1, 82 = 2 and #3 = 3, i.e. the scheme that first appeared in [5], and 
the generalized Fujiwara scheme of order 8 with the choice of parameters 8\ = 1, 82 = 2, #3 = 3 
and 9 4 = 4. 

In order to compare the algorithm to the basic Ninomiya-Victoir scheme we consider an Asian 
call option written on an asset whose price process follows the Heston stochastic volatility model. 
Let X\ be the price process of an asset following the Heston model: 



Xx(t,x) =xi+ / pX 1 (s,x)ds + / X 1 {s,x)yJX 2 (s,t)dB 1 (s) 
Jo Jo 

(21) X 2 (t,x)=x 2 + f a(8- X 2 (s,x))ds 

Jo 
ft 



[3^X 2 (s,t)(pdB 1 (s) + Vl-P 2 dB 2 {s)), 

where x — (x\,x 2 ) S (R>o) 2 , (B 1 (t) , B 2 (t)) is a two-dimensional standard Brownian motion, 
— 1 < p < 1 and a, 8,p are some positive coefficients satisfying 2a8 — 1 > to ensure that the 
volatility does not reach zero. The payoff of the Asian call option on this asset with maturity T 
and strike K is max(Aa(T, x) /T — K, 0), where 

(22) X 3 (t,a) - / X 1 (s,x)ds. 

Jo 

Hence, the price of this option becomes D x E[max(X3(T, x)/T— K, 0)] where D is an appropriate 
discount factor on which we do not focus in this experiment. As in [S] take T = 1, K = 1.05, 
p = 0.05, a = 2.0, [3 = 0.1, 8 = 0.09, p = and x = (1.0, 0.09). 
Up to the error of the magnitude 10~ 6 we have 

£[max(X 3 (T, x)/T - K, 0)] = 6.0473534496 * 10~ 2 

obtained from 8\. Let X(t,x) = {Xi(t, x), X 2 (t, x), X 3 (t, x)) T . SDEs (J21J) and |22]) can be trans- 
formed in the Stratonovich form since X 2 7^ 0: 

2 ft 

X(t,x)=Y ^(A(s,x))od J B l (s), 



i=0 







where 



Vo(yi,y2,V3) = {yi(p~ y - ^-),a(8-y 2 ) - y,2/i) T 



(23) V 1 (y 1 ,y2,y 3 ) = (yiVyi,P0V&>°) T 
V 2 ( yi ,y 2 ,y 3 ) = (0,/3v/(l-p 2 )2/ 2 ,0) T . 

Taking our choice of p = into consideration we get exact solutions of ODEs of the type (|20p 
driven by vector fields V\ and V 2 (see [9] for more details): 

e-xp(tV 1 )(xi,x 2 ,x 3 ) T = (x 1 e t ^,x 2 ,x 3 ), 

(24) / f3t 2 \ 

exp(tV 2 )(xi,x 2 ,x 3 ) T = \ xx,(—+~jxi) ,x 3 J. 
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According to the proof of Theorem [5] we need a Runge-Kutta method of order at least 6 to 
approximate the solution exp(tVo)(xi, X2, x^) T for generalized Fujiwara scheme of order 6 and a 
Runge-Kutta method of order at least 8 for a generalized Fujiwara scheme of order 8 if we want a 
linear algorithm. If we allow quadratic computational cost for the the generalized Fujiwara scheme 
of the weak order 8, it is sufficient to use a Runge-Kutta method of order 4. In our example we 
used 9 stage 7-th order Runge-Kutta method from [11], defined by the Butcher's tableau presented 
in the Appendix. 

The pseudorandom numbers in MC were generated by the Mersenne twister algorithm. The 
QMC was performed using Sobol sequence, generated by the library SobolSeq51.dll provided by 
Broda (see [T]). Both MC and QMC integration were performed using 10 s sample paths. 

The use of exact solutions of ODEs driven by vector fields V\ and Va reduces the computational 
cost of the algorithm by 2Mna YlZ—i 9ki where o designates the order of weak generalized Fujiwara 
scheme divided by 2, M denotes the number of MC/QMC sample paths, n is the number of 
subdivision points and a is the number of operations required for solving ODE's driven by V± or 
V2, if we compare it to the results of Theorem [7] 



Method / n 


2 


3 


4 5 


NV 

GF (order 6) MC 
GF (order 6) QMC 
GF (order 8) MC 
GF (order 8) QMC 


0.00208536744970740 
0.00006154245956983 
0.00005526280089 
0.00004536485526115 
0.0000178413262662 


0.00095536839891733 
0.00003651735446759 

0.0000105789197729 
0.00003694928288030 

0.0000013695959963 


0.00055694952858933 
0.00003522768790512 

0.0000040357269938 0.0000028986604713 
0.000055051968504230 

0.0000010913411477 



Discretization Error vs. Number of Partitions 
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2 3 4 5 6 

Number of Partitions 

Figure 1. Error coming from discretization 

The graph in Fig. [T] clearly shows that the new extrapolation method reduces the order of the 
discretization error in comparison to the original Ninomiya-Victoir algorithm for several magni- 
tudes. In the MC case the discretization error almost immediatly converges to the integration 
error (see Fig [T] and Fig. [2]). Also in the QMC case the discretization error is soon (for small n) 
overshadowed by the integration error caused by QMC integration (see Fig. [J), the weak order of 
the extrapolated algorithms can still be observed from the slope of curves in the graph in Fig. [1] 
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Convergence Error from quasi-Monte Carlo and Monte Carlo 
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Figure 2. Error coming from integration 
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6. Appendix 

We give our original proof of Lemma [2] 

Definition 7. Let g be a Lie algebra. For X, Y E g define c\(X, Y) = X + Y and c n (X, Y) by 
the following recursion formula 

(n + l)c n+1 (X,Y) = ±[X-Y,c n {X,Y)]+ 

+ E A '2 P £ [c kl (X,Y),[...,[c k2p (X,Y),X + Y}...}}, 

p>l,2p<n hi ,. . . ,/c2 p >0, 
feiH \-k 2p =n 

where Kip are coefficients defined in [f4j 2.15.9] 
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For more details about c n (X,Y) see [2J Sec. 2.15]. 
Lemma 4. 

(25) c n (X,Y) = (-l) n+1 c n (Y,X) 

Proof. For n = 1 the assertion is clear. 

Suppose we have c m (X,Y) — (— l) m+1 c m (Y, X) for all m < n. By recursion we obtain 

(n+l)c n+1 {Y,X) 

= ^[Y-X,c n (Y,X)] + 

+ E E [c* 1 (y ) A-),[...,[ Cfcap (y;x),A- + y]...]]. 

fciH Vk-2 p =n 

Using the induction hypothesis and bilinearity of Lie brackets, the above equation transforms into 
(n + l)c n+1 (Y,X) 

= ±(~ir + *lX-Y lCn (X,Y)} + 

+ E K *> E (.-l) kl+ - k2p+2p {c kl (X,Y),{..., 

p> l,2p<n fei ,. . . ,/c2 p >0, 
fclH \-k 2p — n 

[c k2p (X,Y) 7 X + Y]...}} 

= {-lT+^-[X-Y,c n {X,Y)] + 

+ E K *> E K 1 (x,y),[...,[ Cfe2p (x,r),x + y]...]]) 

P>l,2p<M &l,...,fe2p>0, 

fcH hfc2p=" 

= (-ir+ 2 (n+l)c n+1 (X 7 Y) 
which proves the assertion. □ 

Let Ti.d denote ji(q^)- 

Proof of Lemma® The case d = is trivial. Next we consider the case d = 1. Using Baker- 

Campbell-HausdorfT formula to expand tj,i and ji(\og(q^(l))) and applying (|2"5|) proves the for- 
mula Ijlip . By applying Baker-Campbell-Hausdorff formula to the definition of t;^ we get 

n , d = ji(log(qW(d))) =ii(log(exp(log(^l(d-l)))exp(a d ))) 
i i 

= .7; ( E c,£ El 7 ^- 1 ' ad )) ■ 
fe=i j=i 

Suppose that for all n £ N, n < d we have 

, oo 

log(gW(n))=^(-l) i+1 r i)n . 
i=l 

Using Lemma SI the induction hypothesis and the BCH-formula on ji{\og{q [1 \d))) gives us 
ii(l°g(9^(<*))) - j,(log(exp(a d )ex P (log(JW(d- 1))))' 



= jl(Y,c k (A d ,log(qW(d-l)))) 

k=l 
I I 

=i i (E c fc(^.E(- 1 ) i+lr ^-i)) 

fc=l 3=1 
I I 

=j J (E(- 1 ) fe+lc fe(E(- 1 )' +1 - r ^'^))- 
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(26) n(c k C£{~iy +1 r hd -„a d )) = (-l) k+l n (c k (£ Tj , d _ x , a d )) . 

3=1 3=1 



Note that the equality in (|26|) holds trivially for k > I. 

Since Tj^-i is a homogeneous polynomial of degree j, the assertion is clear for k = 1 and all 
Z € N. It is easy to see that for I' < I we have 



(27) j> (cmEl-^V^)) = JV (^(E^ 1 ^. d-i, ad))- 

3=1 3=1 



Let now ^((^.(E^^-l^^d-i,^)) = (-l) fc+ 'ji (c m (X^ = i 7"j,d-i, for all m S {1, . . . ,k} 
and (eN, then we have 



ii(c* + i(x;(-i)'' + S\«i-i>««o) 

3 = 1 



3 = 1 3 = 1 



E ^ £ [ Cfcl (^(-l) 3+1 r^_ liad ), 



P>1 fci,...,fc 2p >0 j'=l 

2p<fc+l fci + ---+fc 2p =/c+l 



3 = 1 

I I 

kEt-^Vi-^Et-^Vi+oi]-"]]) 



3=1 3=1 
1 /l'- 1 ' 



3 = 1 3 = 1 

Z 

+ [ - Od, jj-x(ck(£(-l) i+1 7- A d _!, a d ))] + 

3=1 

i 

+ E ^ E E [^(^(^(-i)^ 1 ^^!,^)), 

P>1 fei,...,/«2p>0 mi,...,m 2 p+l>0 3=1 
2p<fc+l fej^-l [■k 2p =km l J t hm 2 p+i=/ 

I I 

3 = 1 3 = 1 

(-l) m2p+1 + 1 T m2p+1 4-l+3rn 2p+ M] ■■■]])■ 
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Using (|27|) . the induction hypothesis and the bilinearity of Lie brackets the above expression 
transforms into 

-l 



Thus, 



1 /l 

3i{ 

i=i '" ' ~ '" j=i 

i i 

JJ-j( c fc(53 r j,«J-i> a <j))] + - a d ,.7'j_i(cfe(53Tj- i( i_i,a d ))] 

_|_ (— hni2p+l+l+*;iH hfej 

p>l fci,...,fc 2p >0 mi,...,m2p+i>0 
2p<fc+l fcj-l \-k 2p = km 1 -\ i-m 2p+1 =l 

i I 

3 = 1 J=l 
/ 

\_jm 2p { c k 2p i T m 2p +l ,d—l + 3m 2p +i ( a d) 
3 = 1 



Jl 



ii(ck+i(^2(-iy +1 T j: d-i,a d )) 



3 = 1 



I 1 1 
k + l jl ^E^- 1 ~ ad ' Cfc E ^d-i,arf)] 



i i 
+ K 2p Y [ck 1 (Y^ T 3\d-i,ad),[ck*(Y2 l T i,d-i,ad), 

p>l fei,...,fe 2p >0 j=l 3 = 1 

2p<fc+l fc 1 + ...+fe 2 p = fe+l 

i I 

(22 T j,d-li a d),22^j,d-l + a d] ■■ 
3=1 3=1 
I 

= (-l) k+l+X jl{c k+1 (J2r j , d -i,a d )), 

3 = 1 

which is the desired result. □ 
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